
# 读取数据
data <- read.csv("D:\\RLocal\\Class6\\final_merged_data.csv")
# 修改行名为蛋白质名称
rownames(data) <- data$Row.names
data$Row.names <- NULL
# 对数标准化
data2 <- log2(data)
# 提取样本名进行分组
sample_names <- colnames(data2)
group <- gsub("[0-9]+$", "", sample_names)  # 提取字母部分作为分组
group <- as.factor(group)
# 进行ANOVA方差分析
# 设计函数满足题目条件
get_anova_p <- function(protein_expr){
  expr <- as.numeric(protein_expr)
  
  if (length(expr) != length(group) || all(!is.finite(expr))) {
    return(NA)
  } # 如果长度小于理论值或者全为非法值，直接判定为NA
  temp_df <- data.frame(expression = expr, group = group) # 按照分组把数据存入数据框temp_df内
  
  temp_df <- temp_df[is.finite(temp_df$expression), ] # 删除无效数据
  
  if (length(unique(temp_df$group)) < 3) {
    return(NA) # 有效数据少于三个，直接判定为NA
  }
  valid_counts <- tapply(
    temp_df$expression, temp_df$group, length
    )
    
    if (all(valid_counts >= 3)) {
      fit <- aov(expression ~ group, data = temp_df)
      return(summary(fit)[[1]][["Pr(>F)"]][1]) # 返回第一行的p值
    } else {
      return(NA) # 双重保险
    }
  }

# 逐行应用函数求p值 
p_values <- apply(data2, 1, get_anova_p)
# 获得结果
result <- data.frame(Protein = rownames(data2), p_value = p_values)
